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Abstract. A new method, based on the simulated an- 
nealing algorithm and aimed at the inverse problem in 
the analysis of intergalactic (interstellar) complex spec- 
tra of hydrogen and metal lines, is presented. We consider 
the process of line formation in clumpy stochastic me- 
dia accounting for fluctuating velocity and density fields 
(mesoturbulence) . This approach generalizes our previous 
Reverse Monte Carlo and Entropy-Regularized Minimiza- 
tion methods which were applied to velocity fluctuations 
only. The method allows one to estimate, from an observed 
system of spectral lines, both the physical parameters of 
the absorbing gas and appropriate structures of the veloc- 
ity and density distributions along the line of sight. The 
validity of the computational procedure is demonstrated 
using a series of synthetic spectra that emulate the up-to- 
date best quality data. Hi, Cii, Sill, Civ, Siiv, and O vi 
lines, exhibiting complex profiles, were fitted simultane- 
ously. The adopted physical parameters have been recov- 
ered with a sufficiently high accuracy. The results obtained 
encourage the application of the proposed procedure to 
the analysis of real observational data. 

Key words: line: formation - line: profiles - quasars: ab- 
sorption lines - quasars:individual: APM 08279-1-5255 



1. Introduction 

QSO absorption line spectroscopy being a major activity 
at many observatories for the last two decades is now de- 
veloping into a powerful tool for extragalactic research 
thanks to the new generation of large telescopes. The 
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steady improvement in sensitivity and resolution of spec- 
troscopic instrumentation opens new fields in the study of 
QSO absorption systems. It is now becoming possible to 
investigate the intensity fluctuations within the line pro- 
files and thus to estimate hydrodynamic characteristics of 
the absorbing gas. 

The majority of the narrow QSO absorption lines rep- 
resents intervening systems and allows us to probe the 
properties of diffuse matter at very high redshifts. Re- 
solved profiles of hydrogen lines and especially lines of 
heavier elements ('metals') show a diversity of shapes and 
structures. Up to now, their analysis is based on the as- 
sumption that the observed complexity is caused by in- 
dividual 'clouds' randomly distributed along the line of 
sight with slightly different radial velocities. It is also a 
basic assumption that the hydrodynamic ('bulk' or 'tur- 
bulent') velocity distribution inside each cloud is Gaus- 
sian and completely uncorrelated [microturbulence) . This 
model implies that each subcomponent of the complex 
profile being resolved should have a symmetrical profile 
and its shape should not alter with higher spectral res- 
olution. Observations show, however, that the complex- 
ity of the line profiles increases with higher resolution, 
a tendency expected for correlated bulk motions which 
have, in general, non-Gaussian distributions along a given 
line of sight (Levshakov & Kegel 1997; Levshakov, Kegel 
& Mazets 1997; Levshakov, Kegel & Takahara 1999; Pa- 
pers I, II, and HI hereafter, respectively). It follows that 
the microturbulent approximation is not appropriate in 
this case because it does not account for all the relevant 
physical processes involved in the radiative transfer. More- 
over, being applied to real data, the microturbulent analy- 
sis leads to a loss of valuable information contained in the 
observations and may even yield unphysical results (Lev- 
shakov & Kegel 1999; Levshakov, Takahara & Agafonova 
1999; LTA hereafter). The need for more sophisticated 
procedures of data analysis becomes therefore obvious. 

In recent years, it has been shown that accounting for 
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change the interpretation of the Hne measurements sub- 
stantially (Papers I and II). A mcsoturbulcnt approach 
has been already successfully applied to the study of the 
deuterium and hydrogen absorption in Q 1937-1009 (Lev- 
shakov, Kegel & Takahara 1998a), Q 1718+4807 (Lev- 
shakov, Kegel & Takahara 1998b), and Q 1009+2956 
(Levshakov, Tytler & Buries 2000). For all three QSOs 
about the same value for the D/H ratio was derived in 
contrast to the previously announced microturbulent re- 
sults. Our first inversion codes, - the Reverse Monte Carlo 
(Paper III) , and the Entropy- Regularized Minimization 
(LTA), - have been developed for a model of a stochas- 
tic velocity field neglecting any density fluctuations. They 
have been applied to the analysis of the H i and D i lines 
and/or to the metal absorption lines with similar pro- 
files when species trace the same volume elements inde- 
pendently on the density fluctuations. In the present pa- 
per, we extend this study to the inverse problem for a 
model of compressible turbulence when one observes non- 
similar profiles of different atoms and/or ions from the 
same absorption-line system. As in our previous papers, 
we use the term 'turbulence' in a wider sense as compared 
with hydrodynamic turbulence to label the unknown na- 
ture of the line broadening mechanism. In this regard we 
consider any kind of bulk motions (infall, outflows, tidal 
flows etc.) of more or less stochastic nature leading to fluc- 
tuating velocity and density (temperature) fields as con- 
tinuous random functions of the space coordinate along a 
given line of sight within the intervening absorbing region. 

Two noteworthy works have been recently carried out 
aiming at the recovery of the physical intergalactic struc- 
ture from the Lyman-a forest lines. Nusscr & Hachnelt 
(1999a, b) developed an inverse procedure based on the 
relation between density and velocity Fourier coefficients. 
The quality of their recovery is, however, restricted by 
the assumption that the Lyman-a forest structure traces 
mainly the matter density distribution and that the ampli- 
tude of the peculiar velocities is rather small to affect the 
local absorption coefficient significantly. This assumption 
is questionable since there is no simple way to distinguish 
observationally whether the density or the velocity fluctu- 
ations are the main source of the 'line-like' structure ob- 
served in the Lyman-a forest (Levshakov & Kegel 1998). 
Moreover, recent studies of nearby large-scale motions in 
the universe indicate that the Hubble flow is considerably 
perturbed. Peculiar velocities in the range from 300 to 500 
km s""'^ have been found in a sample of galaxies complete 
out to a distance of 60 Mpc (e.g., Watkins 1997; Gramann 
1998; Giovanelli et al. 1998), a fact which should be taken 
into account in the inverse procedures. 

The method described in the present paper is quite 
flexible and equally accounts for the density and velocity 
fluctuations. It is based on a stochastic optimization ap- 
proach similar to that developed in Paper III. We estimate 
simultaneously the physical parameters and appropriate 



butions along the line of sight to reproduce hydrogen and 

metal lines from a given absorption system. In this regard, 
the more spectra of different elements are incorporated in 
the analysis the higher accuracy of the estimation can be 
obtained. 

In §2 our model and the underlying basic assumptions 
are specified. The inversion code is described in §3. The 
validity of the method is tested in §4 using simulated sets 
of noisy line profiles (Hi, Cii, Sill, CiV, Si IV, and Ovi). 
Finally, the main conclusions are outlined in §5. 

2. Hydrogen and metal absorption from QSO 
Lyman-limit systems 

In this section we consider the line formation in a Lyman- 
limit system (LLS), - the intervening absorbing gas being 
optically thin in the Lyman continuum (presumably outer 
regions of a foreground galaxy). To specify the calcula- 
tions, we use the standard photoionization model of Don- 
ahue & Shull (1991, hereafter DS), namely, an extended 
region of thickness L with a given metaUicity. The region 
is ionized by a background photoionizing spectrum given 
by Mathews & Ferland (1987). The gas is assumed to be 
in thermal equilibrium. 

We concentrate our efforts on the LLSs because, on 
one hand, they can be analyzed with a minimum number 
of model assumptions and, on the other hand, they are 
the most promising targets for measuring the extragalac- 
tic deuterium to hydrogen ratio (Buries & Tytler 1998a, b). 
In addition, one may expect that kinematic characteris- 
tics such as dispersions of the velocity and density fluc- 
tuations within LLSs are directly related to the processes 
of galaxy formation. If this is true, these characteristics 
should change with cosmic time (i.e. with z), and we can 
estimate them through the inversion procedure in ques- 
tion. 

The LLSs often show carbon and silicon line absorp- 
tion from different ionization stages and even O VI lines 
(e.g., Kirkman & Tytler 1999). The electron density in 
the LLSs is rather low, rio ^ 10~^ — 10"'^ cm~^, and the 
kinetic temperature T is about 10'* K. Under these condi- 
tions photoionization dominates the ionization structure 
and the relative abundance of different ions of the same el- 
ement is a function of the density only, once the radiation 
field has been specified. 

With the assumed spectral distribution it is convenient 
to describe the thermal and the ionization state of the gas 
in the LLSs through the dimensionless 'ionization param- 
eter' U = nph/riH, equal to the ratio of the number of 
photons with energies above one Rydberg to the total hy- 
drogen density. Following DS, we assume that the ioniz- 
ing background radiation is not time or space dependent. 
Then for a given value of Uph (or the specific radiation 
flux Jo at 1 Ry), the distribution of U along the line of 
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Within the absorbing region the velocity distribution 
v(s) for all species is the same, but their fractional ion- 
izations may vary significantly along the line of sight and 
therefore their apparent profiles may show a diversity of 
shapes. This allows us to tackle the inverse problem, i.e. 
to find such velocity, gas density and temperature distri- 
butions along the line of sight that provide the observed 
variety of the line profiles. 

2.1. Basic equations 

We consider the formation of absorption lines in the light 
of a point-like source of continuum radiation, i.e. the ab- 
sorbing region is assumed to be far from the quasar. The 
observed absorption spectra therefore correspond to only 
one line of sight, i.e. to one realization of the radial veloc- 
ity and gas density distributions. 

For a resolved and optically thin absorption line, one 
observes directly the apparent optical depth r* (A) as func- 
tion of wavelength A 



T*{X)=\n{I,/h) , 



(1) 



where Ix and Jc are the intensities in the line and in the 
continuum, respectively. 

The recorded spectrum is a convolution of the true 
spectrum and the spectrograph point-spread function 0sp 



I 
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(A - A') dA' 



(2) 



where r(A) is the true (intrinsic) optical depth. 

The intensity of the background continuum source Ic 

changes usually very slowly over the width of the spectro- 
graph function, and hence the normalized observed inten- 
sity is 



h 
Ic. 



I 



r(A') , 



,(A- A')dA' . 



(3) 



The instrumental point-spread function (pgp can be de- 
termined experimentally, whereas r(A) is the quantity of 
astrophysical interest. The true optical depth is defined 
through the local absorption coefficient kx{s) by 



t(A) 



/' 



Kx{s)ds = L / Kx{x)dx , 



(4) 



where x = s/L is the normalized coordinate along the line 
of sight. 

In this equation kx{x) is a stochastic variable which 
depends on the random realization of three fields : veloc- 
ity, density and temperature. We write kx as the product 
of the local absorption cross-section per atom, kx{x), and 
the local number density of absorbing atoms, n{x) 



At point X, the absorption cross-section has the form 
kx{x) = ko^x[^Mx),vix)]. (6) 
The quantity ko is a constant for a particular line 



ko 



Tie 
nieC- 



■/absAQ ■ 



(7) 



where me and e are the mass and charge of the electron, 
/abs is the oscillator strength of the line for absorption, 
and Ao is the rest wavelength of the line center. 

The profile function $ describes the local broadening 
at point X. ft is Doppler shifted according to the local 
velocity v{x). Thus we have 

$a[AAd(x), v{x)] = $a{AAd(x), [A - Ao - Ao^]} . (8) 

Under interstellar /intergalactic conditions the profile 
function $ is determined in general by the Doppler effect 
and by radiation damping, i.e. it corresponds to a Voigt 
profile. The Doppler width of the line, A Ad, is determined 
in this case by the thermal width, Vth, 



AAD(a;) = Ao 



vth{x) Ac 



\/2kBT{x)/ms. , 



(9) 



c c 

and hence vth characterizes the width of the local absorp- 
tion coefficient at a given value of v (here /cb is Boltz- 
mann's constant, ma is the mass of the ion under consid- 
eration, and T{x) is the local kinetic temperature). 

Now, let n^ix) be the local number density of element 
'a', na,i(a;) be the density of ions in the ith ionization 
stage, and Za = ria/^H be the relative abundance of this 
element which is regarded to be a constant along the line 
of sight. According to the standard model of DS, the frac- 
tional ionization of ions {a, i} 

ne,,iix) 



T — 



na.{x) 



(10) 



may be described by a function of U only, Taj = 

Then for a fixed value of A within the line profile, the 
optical depth is given according to equations (4)-(6) by 

Ta,i(A) = 

koL [ n^{x)T^4U{x)]<S>x[AXj,{x),v{x)]dx, (11) 
Jo 

or with na,{x) = Za,nu{x) 

Ta,i(A) = 

koZ^L [ nHix)r^4U{x)]^x[AXD{x),vix)]dx. (12) 
Jo 

Having defined nii{x), the ionization fractions Ta,i can 
be calculated using a photoionization model correspond- 
ing to that of DS or Ferland (1995). In the present paper 
we used the results of DS since we may neglect the ra- 
diation transfer and high-density effects incorporated in 
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Fig. 1. Velocity plots of Si iv (a) and Si ill (b) synthetic profiles 
calculated for the = 3.514 system towards QSO 08279-1-5255 
by Levshakov, Agafonova & Kegel (2000). The underlying 
velocity, density and ionization parameter distributions are 
shown in panels (c), (d), and (o), respectively. The diversity 
of fractional ionizations in the marked areas is illustrated in 
panel (f) 



(see e.g. Tytler et al. 1996). This fitting method would 

yield physically reasonable results for gas in spatially sep- 
arate homogeneous clouds if we could know a priori the 
actual number of these clouds and have a justification for 
their homogeneity. As far as the number of subcomponents 
is unknown, the VPF method loses uniqueness and stabil- 
ity. The main question remains whether the complex line 
shape is caused by separate homogeneous clouds or by a 
continuous medium with fluctuating characteristics. There 
are no unambiguous observational data favoring one or the 
other assumption. But some indirect information such as 
the increasing complexity of the line profiles with increas- 
ing spectral resolution and recently performed 3D hydro- 
dynamic calculations (see e.g. Ranch, Haehnelt & Stein- 
metz 1997) give higher weight to the model of a fluctuat- 
ing continuous medium. Applied to the inversion problem, 
this has important consequences. 

Fig. 1 illustrates step-by-step the line formation pro- 
cess in a stochastic medium. The synthetic spectra of 
Si IVA1393 (panel a) and Si iiiA1206 (panel b) are formed 
in an absorbing region with the fluctuating velocity and 
density fields shown in panels c and d, respectively. The 
corresponding distribution of the ionization parameter 
U {x) is shown in panel e. 

Let us consider a point within the Si iv and Si iii pro- 
files, for example at the radial velocity v = Q km (pan- 
els a and b). Then filled circles in panel c mark those 
volume elements having the same radial velocity. The cor- 
responding points in panel d show the density in the cor- 
responding areas contributing to the 'observed' intensity. 

It is clearly seen that for a given point within the 
line profile the observed intensity results from a mixture 
of different ionization states (labeled portions of the U- 
distribution in panel e). It is also obvious that each point 
of the line profile is caused by a different combination of 
these states. As an example, panel f shows the different 
fractional ionizations of Si III and Si IV at the points with 
w = km s~^. 

In practical applications, the metal abundances, Zg_, 
are usually estimated from the ratio of the total column 
densities A^a and A^h by means of the 'correction for the 
ionization' : 



iVa,i THi(t/) 



(13) 



2.2. Some aspects of line formation in clumpy stochastic 
media 

When using the conventional Voigt-profile fitting (VPF) 
procedure several subcomponents are usually assumed to 
describe an absorption line which is too complex to be 
fitted by a single Voigt profile. Each component is sup- 
posed to represent a homogeneous plane-parallel slab of 
gas with its own physical parameters : radial velocity, the 



The mean value of U is usually estimated from the 
ratio Ng,,i/Ns,j. While this procedure is suggestive, it is 
incorrect if the degree of ionization varies along the line 
of sight. From eq. (12) and accounting for J <1?a dA = 1, it 
follows that 



J nnix) Ta, j [U{x)] dx f^.j 



(14) 



where Ta,i and Ta,j denote the mean density-weighted 
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The ratio Ta,i/Ta,j as a function of U can be calcu- 
lated from the adopted photoionization model. Then the 
mean value of the ionization parameter, [/, is estimated 
using the implicitly assumed relation 



(15) 



Since Ta,j is a non- linear function of uh, it is clear 
that the equality T = T(C7) holds only in the case of a 
imiform density. Applying relation (15) to the case of a 
fluctuating gas density leads to a wrong value oiU - and 
moreover the estimation of the mean ionization parameter 
will depend on the specific pair of ions chosen to calculate 
the left-hand side of (14). 

To illustrate the numerical differences that may oc- 
cur between the micro- and mesoturbulcnt approaches 
we consider the following example. The absorption lines 
shown in Fig. 1 have approximately the same total col- 
umn densities [the data are taken from Ellison et al. 
1999 (logA^siiv - 12.8 cm-2) and Molaro et al. 1999 
(logATsini ~ 12.8 cm-2)]. The ratio TVsini/^siiv 1 sug- 
gests that log U ~ —2.2, i.e. ?7 ~ 6.3 x lO^"*. However, the 
mesoturbulent solution yields a 2.5 times higher value for 
the mean ionization parameter, i.e. C7 ~ 1.6 x 10~^ (for 
more details, sec Levshakov, Agafonova & Kegel 2000). 

It is worthwhile to conclude this section with a few re- 
marks concerning the VPF method. Prom a mathematical 
point of view the VPF is quite consistent since it simply 
means the representation of the unknown radial velocity 
distribution function by the sum of several Voigt func- 
tions. In such a way we can very well reproduce the line 
profiles and, in the case of unsaturated lines, calculate ac- 
curately the total column densities. But if we wish to go 
further and estimate physical parameters like kinetic tem- 
perature, ionization and metallicity, the VPF may pro- 
duce misleading results. If reality corresponds to a model 
like that shown in Fig. 1, it is physically unjustified to 
interpret each subcomponent as belonging to one distin- 
guished cloud since as shown in Fig. 1 the contribution 
to any point within the line profile does not come from a 
single separate area but from volume elements distributed 
over the whole absorbing region. 

3. The Monte Carlo Inversion 

The Monte Carlo Inversion (MCI) of line profiles includes 
the estimation of the physical parameters and simultane- 
ously of three random functions nH(a;), T{x), and v{x). For 
the equilibrium models of DS, the kinetic temperature is 
determined from the energy balance equation in the opti- 
cally thin limit. In this case T depends on the gas density 
only, the dependence being described approximately by a 
power law, i.e. 



For an individual realization the remaining two random 

fimctions n-aix) and v{x) are determined as follows. 

The hydrodynamic velocity and density fields are 
formed as a result of interference of many independent and 
random factors infall and outflows, rotation, tidal flows, 
shock waves etc. In this case, we may consider the fluc- 
tuating amplitude of the velocity or the density along the 
line of sight as a kind of Brownian motion which is math- 
ematically described by the so-called Markovian processes 
(e.g. Rytov et al. 1989). In particular, in case of Gaussian 
fields the finite difference representation of the velocity as 
a Markovian process has the form 



v{x + Ax) = /v v{x) + e{x + Ax) 



(17) 



where /v = Rv{Ax)/al, Ry being the correlation between 
the velocity values at points separated by a distance Ax, 
i.e. i?v = {v{x + Ax)v{x)), the velocity dispersion of 
bulk motions, and e a random normal variable with zero 
mean and dispersion 



(18) 



The density and velocity are in general related through 
the continuity equation. However, since we consider one 
component of the velocity field only, the correlation be- 
tween this component and the density is less tight and 
in lowest order approximation we may consider the two 
quantities to be statistically independent of each other. 

Suppose further that the logarithmic density distribu- 
tion is normal and that the dimensionless variable y{x) = 
nii{x)/no, with no being the mean gas density has the 
second central moment CTy characterizing the dispersion 
of the density field along the line of sight [this definition 
implies that the expectation value E{y) = 1 for any x]. 
Then consider, by analogy with v{x), a random function 
u{x) with the following representation 



^{x + Ax) — fi, v{x) + e{x + Ax) , 



(19) 



where = Ru{Ax)/a^, R„ being the correlation between 
the values of u at points separated by a distance Ax, 
the logarithmic density dispersion, and e a random normal 
variable with zero mean and dispersion 



(20) 



Having defined y{x), the distribution of the hydrogen 
density can be obtained through Johnson's transformation 
(Kendall & Stuart 1963) : 

nH(ar)=noe^(-)-^-^ (21) 
with the relationship between the dispersions and Uy 
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Fig. 2. Model A, Markovian fields : the normalized velocity (a), density (b), kinetic temperature (c), and ionization parameter 
(d) distributions vs space coordinate (T is given in units of Iff* K). The corresponding density- weighted distribution functions 
for the radial velocities of the gas and the individual ions are shown in panels (e) and (f). Panels (g) -(1) represent simulated 
spectra convolved with an instrumental profile of FWHM = 7 km s"'^ and added Gaussian noise with S/N = 50 (dots with 
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With (21), we can rewrite the expression (12) for 
Ta,i (A) in the form 



X / e'^(^)Ta,i[C/(x)]$A[AAD(x),i;(x)]dx. 
Jo 



(23) 



where A^'o = no L is the expectation value of the total 
hydrogen column density. 

Equation (23) can be simplified if we have a fully re- 
solved profile of an optically thin line. In this case the total 
ion column density along the given line of sight, N^^^i, can 
be found directly from the observed profile (Paper II), and 
thus we can eliminate the unknown abundance from (23). 
Integrating (23) over A yields 



i.o 
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(24) 



where the mean density-weighted fractional ionization 
Ta,* is 



Ta,i = e-^-' e''(-)Ta,i[f/(a;)]da;. 
Jo 

Inserting (24) into (23) leads to 



„.(A) = 



(25) 



(26) 



Fig. 3. Model B : adopted correlation function of the velocity 
fluctuations (dotted line) and corresponding correlation func- 
tion of the density contrast fluctuations (solid line) 



column density), the mean value of the ionization param- 
eter Uo, the rms values of the density and velocity fluctu- 
ations, cTy and av, respectively, the correlation coefficients 
fi, and /v, and the abundances Zg, of the elements involved 
in the analysis. All these parameters are components of 
the parameter vector 9. Additionally to these physical pa- 
rameters, the density and velocity distributions are to be 
known. The continuous random functions v{x) and ^{x) 



i<T^ / i/(x) -v rrr/ M;iirA\ ^'^ /mj are represented in our computations by their values ii,- and 
2'^. e ^^^Ta,i ?7(a;) $A AAD(a;),u(a;) da;. , j . .i, n j i ^ ■ . mi, 

Jo Uj sampled at the equally spaced spatial pomts Xj. They 



Next, the obvious relationship 



U{x) = 



nii{x) no 



(27) 



are computed with the relations (17) and (19) which cor- 
relate the values at neighboring points (j = 1, 2, k). To 
estimate 9, {vj}^^^ and {yj}^^i from the absorption line 
profiles, we have to minimize the objective function 



shows that the mean ionization parameter Uq is given by 
the expression 



-I K Pf: 

^ e=i i=i 



■obs 



(29) 



{U{x)) =Uo^ ^e^' 
no 



(28) 



(use again the Johnson transformation to evaluate the 
probability density distribution ofU). Here M is a multi- 
plicative factor introduced by DS in oder to express nph in 
terms of Jo- For the photoionizing spectrum of Mathews 
& Ferland, M is 5.19 x 10^^ (in cgs units). 

The result (28) has quite interesting consequences. 
Namely, if the density is completely homogeneous (cry = 0) 
then Uo = MJo/no- But if the density field is perturbed 
{ai, > 0), then with the same mean density no and the 
same background radiation flux Jo one obtains a higher 
value of Uq without any additional sources of ionization. 
Intermittent regions of low and high ionization caused by 
the density fluctuations may occur in this case along a 
given line of sight. 

Now we can describe our model in detail. It is fully de- 



Here J^£^Xi is the observed normalized intensity of spectral 
line £ according to eq. (3), and ag^i is the experimental 
error within the ith pixel of the line profile. -F^^' (9) is the 
simulated intensity of line £ at the same ith pixel having 
wavelength Aj. The total number of spectral lines involved 
in the optimization procedure is labeled by and the 
total number of data points P = X^^Li ^ii where Pi is 
the number of data points for the ^th line. The number of 
degrees of freedom is labeled by p. 

Let us go on to consider the correlation coefficients /v 
and fv As known, the Markovian processes have exponen- 
tial correlation fimctions : R{Ax) oc cxp(— A.x/Z), where 
I is referred to as the correlation length and depends, in 
general, on the spatial scale. For the step size Ax /I <C 1, 
the correlation coefficient / is very close to unity. The 
computational procedure is quite insensitive to the con- 
crete values of / because of the stochastic relations (17) 
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the same manner as other physical parameters. It appears 
to be more suitable to fix a few sets of /v and fi, and then 
to carry out the estimation of other parameters with each 
of these sets. Hence, in reality the parameter vector 6 does 
not include components /v and f,^. 

The computational scheme is similar to that described 
in Paper III, but in the present work we used instead of the 
classical Metropolis method the more efficient generalized 
simulated annealing algorithm (Xiang et al. 1997). 

Model calculations with different synthetic spectra 
have shown that the minimization of in form of eq. (29) 
does not allow us to recover the physical parameters with 
a sufficiently high accuracy (see also LTA). To stabilize 
the estimation of 9, we need to augment by a reg- 
ularization term (or a penalty function). The choice of 
the penalty function is rather heuristic and depends on 
the particular problem. For instance, when we try to fit 
model spectra to real noisy data it is reasonable to stop 
the minimization procedure at the value of — 1 in order 
to avoid fitting intensity fluctuations caused by the noise. 
Another restriction stems from the restoring procedure for 
the density and velocity distributions. If we use instead of 
the processes (17) and (19) their normalized analogies : 
v*{x) = v{x)/ay and i'*{x) = i'{x)/ai, with the disper- 
sions fTv* = <Jv* = 1, more stable estimations for and 
(T,y arc obtained. In addition, the process //* should be 
centered, i.e. (v*) =0 (the center of the velocity distri- 
bution is determined by the relative positions of spectral 
lines and in principle can differ from zero) . Accounting for 
all these remarks, the modified objective function is now 
written in the form 

£ = Ix' - 1| + a(|(^*)| + - 1| + |av. - 1|) , (30) 

where a is a constant of order unity and i)* , ct^., a^' are 
the quantities derived from the current configuration. It 
follows that in the vicinity of the global minimum we must 

have £min ~ 0. 

If the optimization of C gives {v*) ^ 0, the estimated 
rms velocity ay is biased. To correct its value, the following 
relation can be used 

a^' = CTy (^l-k{vr /p<.v;r^ . (31) 

The practical implementation of the proposed compu- 
tational procedure for recovering the physical parameters 
from spectral lines of different ions is described in the fol- 
lowing section. 

4. Numerical Experiments 

We tested the numerical procedure described above by in- 
verting synthetic spectra with known physical parameters 



Table 1. Model parameters used to generate synthetic spectra 



model 


no 


Jo 


L 








cm~* 


ergs/s cm^ Hz 


kpc 


km/s 




A 


0.001 


1.510"^^ 


5.5 


25 


1.5 


B 


0.005 


3.0 10"^^ 


1.5 


20 


0.8 



the line of sight. We have considered two models A and B 
with the physical parameters listed in Table 1. 

To calculate random fields several procedures can 
be utilized. The most simple one has been realized in 
model A : the velocity and gas density fields were pro- 
duced by means of two independent Markovian processes 
according to eqs. (17) and (19) with /v = 0.995, = 
0.999 and Ax = 1/300. The fractional ionization for the 
different ions included in the analysis were then calcu- 
lated using the model of DS. It is worth emphasizing that 
this model may not necessarily capture the physics most 
correctly since it assumes a specific type of the ionizing 
radiation, but since our purpose in this paper is to show 
the principal possibility of recovering the adopted physi- 
cal characteristics of the absorbing region from the syn- 
thetic spectra, the DS model was chosen because of its 
transparency and easiness to compute. The synthetic spec- 
tra were then calculated according to eq. (12) with fixed 
metallicity for all ions Z = 0.1 Zq (solar abundances were 
taken from Grcvesse 1984). To complete the procedure 
the spectra were convolved with a Gaussian point-spread 
function having FWHM = 7 km s^^ and a Gaussian noise 
with S/N = 50 was added to each pixel (dots with error 
bars in Fig. 2). 

Fig. 2 shows the hydrodynamic fields (normalized ve- 
locity u = v/(7y and gas density) together with the kinetic 
temperature (in units of 10^ K) and the ionization pa- 
rameter U/Uo distributions along the line of sight (panels 
a - d) . The histograms in panels e and f give the density- 
weighted velocity distributions of the different ions which 
determine the shapes of the spectral lines presented in 
turn in panels g - 1 for the H i Lyl2, C il, C IV, Si ii. Si iv, 
and O VI lines, respectively. The Lyl2 line was chosen for 
the analysis because it is unsaturated and, hence, allows 
us to estimate the neutral hydrogen column density quite 
accurately. The choice of metallic ions was motivated by 
the fact that their lines are formed in different parts of 
the adsorbing region : Cii, for example, needs denser gas 
and lower temperature to be produced (the area in the 
vicinity of a; ~ 0.2), whereas Ovi mainly is found in the 
very diffuse and hot areas {x ~ 0.02, 0.38 and 0.8). Such a 
diversity of lines in question enables us to restore the un- 
derlying density and radial velocity distributions in more 
detail. 

Model B differs from model A not only by the values of 
physical parameters but by the method to generate ran- 
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Fig. 4. As Fig. 2, but for model B, non-Markovian fields 
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density contrast, 6{x), distributions, the corresponding 

power spectra were utihzcd. The computational scheme 
was similar to that proposed by Bi, Ge & Fang (1995) 
but with another set of power spectra. We assume that 
the ID power spectrum for the velocity is given and then 
we calculate the ID power spectrum for 6. In particular, 
for the normalized radial velocity distribution, u{x), the 
following ID power spectrum was chosen [it corresponds 
to the correlation function Ci^e{x) adopted in LTA] : 



-Ik 



2 r(2 -£)/<= 



-2 ' 



(32) 



where T is the gamma function, I and e are fixed constants 
(e < 2 and I is in units of L). 

We can calculate the ID power spectrum for the den- 
sity contrast from its 3D spectrum which in turn is defined 
by the 3D power spectrum for the velocity : Sg^ = S^^ . 
(This relation stems from the continuity equation when it 
is linearized in comoving coordinates). The 3D power spec- 
trum is related to its ID spectrum via (see Monin & 
Yaglom 1975) : 



27rfc2 
and thus 



dfc' 



27rfc2 



2 5'i"(fc)-fc 



dfc 



(33) 



(34) 



The ID power spectrum for the density contrast is then 
determined directly from its 3D spectrum (Monin & Ya- 
glom) : 



5f (fc')fc'dfc' 



(35) 



Given the power spectra S'^^(fc) and Sg^{k), which 
are the Fourier transforms of the covariance functions of 
u and S, both random fields u{x) and S{x) can be modeled 
using the moving-average method as described in LTA. In 
order to avoid the appearance of negative n{x), the gas 
density itself is assumed to have a log-normal distribution 
(see e.g. Bi, Borner & Chu 1992; Bi & Davidsen 1997). It 
is calculated from S according to 



n{x)/no = e' 



(36) 



The density contrast, as calculated from (36), equals 5{x) 
in linear approximation. 

The correlation functions of the processes u{x) and 
5{x) with the foregoing power spectra and with the 
adopted parameters e = 0.1 and / = 0.1 are shown in 
Fig. 3. They differ significantly from the exponential form 
corresponding to the pure Markovian process. So the ob- 
jective to utilize these kind of stochastic fields was to test 
whether we can approximate a finite realization of an ar- 



corresponding hydrodynamic fields and line profiles for 
model B are shown in Fig. 4 (notations are the same as 
in Fig. 2). 

The recovering procedure aimed at the minimization of 

the objective fimction C, eq. (30), consists of two repeating 
steps : firstly a random point in the physical parameter 
box is chosen and then the appropriate optimal velocity 
and density configurations arc estimated. (In our actual 
calculations we used instead of the physical parameters 
A^o and Uq their reduced values Nq = Nq e~ a'^" and Uq = 
?7oc-2<^,.. This was made to achieve better stability of the 
results) . 

In both steps the stochastic acceptance rule was em- 
ployed. The move in the parameter space or in the 
velocity-density configuration space is accepted if it leads 
to a lower current value of C. If C jumps up, then the move 
might be accepted according to an acceptance probability 
Pa- Two types of the acceptance probability were used in 
the procedure : the standard Boltzmann-type statistics as 
utilized in the Metropolis algorithm 



Pa 



-AC/T 



(37) 



and the Tsallis statistics being the generalization of the 
Boltzmann distribution (e.g. Xiang et al. 1997) 



Pa = [1 - (1 - ga)A£/r] 



l/(l-9a) 



(38) 



where q^, is a fixed constant and T is the so-called anneal- 
ing parameter ('artificial temperature'). 

In general the acceptance after Tsallis leads to quicker 
convergency, but it depends strongly on the value of (^a 
and requires a few initial trial computations to adjust 
this parameter. The annealing parameter T should de- 
crease with the simulation time to ensiire the convergence 
to the global minimum. We used the 'cooling' rate of 
type T = 7o/(l -|- Ngteps)'', with the initial value % and 
the power index q having different values for both proce- 
dure steps (i.e. the move in the parameter space and the 
move in the density- velocity configuration space). Their 
values were adjusted in numerous numerical experiments 
with the aim to reach better stability and quicker con- 
vergency. The results presented below were obtained with 
TqP'^'' = 2.0, = In(x^), <?par = 0.5, gd-v = 2 and 

ga = — 1 (here Xo is the value for the initial step in 
the parameter space). The visiting distribution (a local 
search distribution providing trial values) was Gaussian 
in all cases with the fixed dispersion of 0.1 (in dimension- 
less units). 

The results of the computations are shown in Fig. 2 
(model A, Markovian fields) and in Fig. 4 (model B, non- 
Markovian fields) by solid lines (panels g-1). For both 
cases the step size was Aa; = 1/100 and the correlation 
coefficients fv = = 0.95. These values were chosen as 
the optimal ones after several trial runs with different Aa;, 
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Fig. 5. MCI reconstruction of the velocity and density fields for three different runs (al,a2), (bl,b2), and (cl,c2) for model A, 
and (el,e2), (fl,f2), and (gl,g2) for model B (different types of dashed and dotted lines). The adopted distributions are shown by 

solid lines. Panels (d) and (h) represent the corresponding density-weighted velocity distribution functions (solid line histograms 
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Table 2. Physical parameters derived by the MCI method (rows 1-11) from the synthetic spectra shown in Fig. 4. The first 
and the last row lists, respectively, the adopted values (model B) and the median estimations 



Uo No (Tv (Ty Zc Zsi Zo 5"^* 

cm~^ km/s 





4, 


.426e- 


■3 


1 


.583el9 


20.0 


0.8 


4. 


.900e-5 


3, 


.550e- 


■6 


8. 


.lOOe-5 












1 


4 


.392c- 


■3 


1 


.819cl9 


20.04 


0.722 


4, 


.839C-5 


3 


.488c- 


-6 


5, 


.318C-5 


1 


.057 


-0.008 


1.020 


1.025 


2 


3, 


.871c- 


■3 


1 


.744cl9 


22.95 


0.729 


4, 


,835c-5 


3, 


.558c- 


-6 


8, 


.409c-5 


1 


.031 


0.033 


1.017 


1.006 


3 


5, 


.532c- 


■3 


1 


.580el9 


24.57 


1.215 


4, 


.853C-5 


3 


.540c- 


-6 


4, 


,516c-5 


1 


.046 


-0.028 


1.009 


0.937 


4 


3 


.878e- 


■3 


1 


.727el9 


21.08 


0.803 


4, 


.974e-5 


3 


.602e- 


■6 


8, 


.336e-5 


1 


.095 


4.0e-4 


0.993 


1.009 


5 


4, 


.081e- 


■3 


1 


.837el9 


21.97 


0.695 


4, 


.936e-5 


3 


.490e- 


■6 


7, 


.794e-5 


1 


.100 


-0.013 


1.005 


0.986 


6 


3 


.104c- 


■3 


1 


.964cl9 


22.92 


0.447 


4, 


.919e-5 


3 


.456c- 


-6 


1, 


,406c-4 


1 


.095 


-2.0C-4 


1.012 


0.966 


7 


5, 


.683e- 


-3 


1 


.490el9 


20.98 


1.419 


4, 


.865e-5 


3 


.628e- 


-6 


5, 


.097e-5 


1 


.103 


-0.019 


0.991 


1.001 


8 


4, 


.083e- 


■3 


1 


.780el9 


19.56 


1.206 


4. 


.988e-5 


3, 


.458e- 


■6 


7. 


.769e-5 


1 


.093 


-0.068 


0.983 


1.096 


9 


4, 


.520e- 


■3 


1 


.546el9 


22.77 


1.012 


4. 


.840e-5 


3, 


.701e- 


■6 


7. 


.674e-5 


1 


.151 


-0.004 


1.002 


0.991 


10 


4, 


.980e- 


■3 


1 


.630el9 


26.05 


0.998 


4. 


.913e-5 


3, 


.516e- 


■6 


7. 


.615e-5 


1 


.155 


0.008 


1.014 


0.994 


11 


3, 


.848e- 


■3 


1 


.821el9 


20.28 


0.679 


4. 


.961e-5 


3, 


.535e- 


■6 


8. 


.471e-5 


1 


.090 


-0.013 


1.000 


1.017 




4, 


.083e- 


•3 


1 


.727el9 


21.97 


0.803 


4. 


.913e-5 


3, 


.535e- 


■6 


7. 


.674e-5 













The MCI is a probabilistic procedure. It follows that 
one does not know in advance if the global minimum of 
the objective function is reached in a single run. Several 
runs with different initial values should be performed to 
obtain the parameter estimations and their confidence lev- 
els. To illustrate the probabilistic character of the MCI 
solutions, Fig. 5 shows the results of three different runs 
for model A (left-side panels) and for model B (right-side 
panels). Solid lines in all panels represent adopted fields, 
whereas dashed and two types dotted lines are the MCI so- 
lutions. All recovered fields yield ^ 1) implying that the 
absorption spectra calculated in different runs are indis- 
tinguishable from each other. The example shown by solid 
lines in Fig. 2 or 4 (panels g - 1) is just one of the solutions 
obtained for the corresponding model A or B. However, it 
is not possible to restore the exact pixel-to-pixel structure 
of the stochastic velocity and density patterns along the 
line of sight - many configurations are possible with equal 
probability. But we see that all these configurations have 
the same (statistically indistinguishable) density-weighted 
velocity distributions (panels d and h), since as already 
mentioned above only these distributions determine the 
observed line shapes. It is obvious that the quality of re- 
covering depends crucially on the diversity of ion profiles 
accounted for in the MCI analysis : the more ions origi- 
nating from regions with different physical conditions are 
considered the higher reliability of the estimated parame- 
ters can be reached. 

Table 2 lists the results of 11 nms for model B 
(model A shows entirely analogous results). It is seen that 
the accuracy of the estimated parameters is rather high : 
the median values do not differ from the adopted ones by 
more than 10%. As for the confidence levels, it is obvious 
that much more runs are needed to calculate them with a 



tations since in the present context the exact values of the 
confidence levels are not of importance. But from the runs 
performed we can make some qualitative conclusions. For 
instance, the metallicities are estimated with very high ac- 
curacy (the error in Zq is explained by the weakness of the 
oxygen line), i.e. the numerical procedure is quite sensitive 
to these parameters. The dispersions of and cTy show 
significantly wider margins caused by the randomness of 
the fields which they characterize. The estimations of Uo 
and A^o have an accuracy of about 20% that can con- 
sidered as quite acceptable for these physical parameters. 
A higher accuracy can be obtained for the estimations of 
the total column densities of the individual ions and for 
the mean (density-weighted) kinetic temperatures 



J ni,^i{x) dx 



(39) 



These physical parameters are listed in Table 3 in the same 
order as the results presented in Table 2. 

To complete this section, we note that the computa- 
tional time of one run with the adopted set of the spec- 
tral lines and with two 100 components vectors represent- 
ing the velocity and density random configurations takes 
about 10 hours on a Pentium III processor with 500 MHz. 

5. Conclusions 

We have developed a new method to solve the inverse 
problem in the analysis of intergalactic (interstellar) hy- 
drogen and metal lines arising from clumpy stochastic 
media. In the method, the random velocity and density 
configurations along the line of sight are approximated 
by Markovian processes. The global optimization method 
based upon simulated annealing is then used to fit theo- 



S.A. Levshakov et al.: Monte Carlo inversion of line profiles 



13 



Table 3. Physical parameters derived by the MCI method (rows 1-11) from the synthetic spectra shown in Fig. 4. The first 
and the last row lists, respectively, the adopted values (model B) and the median estimations 



Hi 


Cii Civ Sin Siiv 


Ovi 


Hi Cn Civ Sin Siiv 


Ovi 




Column densities in cm~^ 




Mean kinetic temperatures in 10* 


K 





5.92el6 


5.52el3 


2, 


.76el4 


4.43el2 


1 


.80el3 


6.22el2 


1 


.75 


1, 


.74 


1, 


.88 


1, 


.74 


1 


.84 


2. 


.16 


1 


5.95el6 


5.24el3 


2 


.86el4 


4.33el2 


1 


.77el3 


5.14el2 


1 


.74 


1, 


.72 


1 


.93 


1, 


.73 


1 


.87 


2. 


.14 


2 


6.07el6 


5.61el3 


2 


.65el4 


4.55el2 


1 


.75el3 


5.90el2 


1 


.74 


1, 


.72 


1, 


.89 


1, 


.73 


1 


.84 


2. 


.14 


3 


5.96el6 


5.52el3 


2 


.77el4 


4.45el2 


1 


.77el3 


6.05el2 


1 


.75 


1, 


.73 


1, 


.91 


1, 


.73 


1 


.85 


2, 


.32 


4 


5.85cl6 


5.50ol3 


2 


.79cl4 


4.40el2 


1 


.79el3 


5.55ol2 


1 


.74 


1, 


.72 


1, 


.89 


1, 


.73 


1 


.85 


2, 


,14 


5 


6.06el6 


5.68el3 


2 


.86el4 


4.43el2 


1, 


.76el3 


6.64el2 


1, 


.74 


1, 


.72 


1 


.91 


1, 


.72 


1, 


.86 


2. 


.12 


6 


5.98el6 


5.56el3 


2 


.82el4 


4.33el2 


1 


.78el3 


6.13el2 


1 


.74 


1, 


.73 


1 


.87 


1, 


.74 


1 


.84 


1. 


.97 


7 


5.96el6 


5.51el3 


2 


.74el4 


4.53el2 


1 


.78el3 


6.30el2 


1 


.73 


1, 


.72 


1 


.92 


1, 


.72 


1 


.86 


2, 


,31 


8 


5.96el6 


5.65el3 


2 


.87cl4 


4.34el2 


1 


.77el3 


6.07el2 


1 


.75 


1, 


.74 


1, 


.89 


1, 


.74 


1 


.85 


2, 


,21 


9 


5.91el6 


5.47el3 


2 


.55cl4 


4.61el2 


1 


.75el3 


6.48el2 


1 


.73 


1, 


.72 


1, 


.89 


1, 


.73 


1 


.84 


2, 


,24 


10 


5.97el6 


5.58el3 


2 


.79el4 


4.42el3 


1, 


.76el3 


8.73el2 


1, 


.74 


1, 


.73 


1 


.90 


1, 


.73 


1, 


.85 


2. 


,26 


11 


5.88el6 


5.55el3 


2 


.83el4 


4.37el3 


1 


.78el3 


5.98el2 


1 


.75 


1, 


.75 


1, 


.90 


1, 


.74 


1 


.85 


2. 


,10 




5.96el6 


5.55el3 


2 


.79el4 


4.42el3 


1 


.77el3 


6.07el2 


1 


.74 


1, 


.72 


1, 


.90 


1, 


.73 


1 


.85 


2. 


,14 



The proposed procedure allows us to estimate the 

physical parameters of the absorbing gas such as column 
densities, metal abundances, mean (density-weighted) ki- 
netic temperatures for each ion, and mean ionization pa- 
rameter together with the hydrodynamic characteristics 
- the radial velocity dispersion and the dispersion of the 
density fluctuations. 

The computational scheme has been tested on a variety 
of synthetic spectra that emulate modern observational 
data : the absorption lines of Hi, Cii, Siii, Civ, Siiv, 
and O VI which arc usually observed in the Lyman limit 
systems (iVni < 3 x 10^'' cm~^). The ionization structure 
of the absorbing region was calculated using the standard 
photoionization model of Donahue & ShuU (1991) with a 
background ionizing spectrum given by Mathews & Fer- 
land (1987). 

The inversion procedure proved to be very effective and 
robust allowing us to recover the physical parameters with 
reasonable accuracy albeit the structure of the random ve- 
locity and density fields cannot be restored with a pixcl- 
to-pixel conformity. However, the integral characteristics 
of these random fields, namely, the density- weighted veloc- 
ity distribution, can be estimated quite precisely. Thus we 
can conclude that our procedure provides reliable results 
and can be applied to the analysis of real data. 

Note that while performing the inversion of absorp- 
tion lines, one has to take into account the following. All 
our computational tests have been carried out under the 
assumption that the spectrum of the ionizing radiation 
is known, i.e. we used the same Mathews & Ferland spec- 
trum to generate 'observational' data and to fit them with 
our theoretical profiles. In reality, the characteristics of 
the ionizing radiation are not known exactly. Therefore 
in real applications several types of the background pho- 



computational results are aff'ected by different types of 
the background ionizing radiation will be studied in detail 
elsewhere. 

The proposed method has been successfully applied to 
the analysis of QSO high resolution spectral data with 
possible deuterium absorption at = 3.514 towards 
APM 08279-F5255 (Levshakov, Agafonova & Kegel 2000). 
It has been demonstrated that the blue-side asymmetry 
of the hydrogen Lya line can be explained quite naturally 
by an asymmetric configuration of the velocity field only. 
The results obtained revealed a considerably lower neutral 
hydrogen column density as compared with the VPF mea- 
surements performed by Molaro et al. (1999). In contrast 
to Molaro et al., we have managed to fit simultaneously all 
absorption lines observed in this system. These results can 
be considered as encouraging and favor the application of 
the developed computational procedure to the analysis of 
other high quality observational data. 
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